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Abstract. The ubiquity of oscillations in epidemics presents a long standing 
challenge for the formulation of epidemic models. Whether they are external and 
seasonally driven, or arise from the intrinsic dynamics is an open problem. It is known 
that fixed time delays destabilize the steady state solution of the standard SIRS model, 
giving rise to stable oscillations for certain parameters values. In this contribution, 
starting from the classical SIRS model, we make a general treatment of the recovery 
and loss of immunity terms. We present oscillation diagrams (amplitude and period) 
in terms of the parameters of the model, showing how oscillations can be destabilized 
by the shape of the distributions of the two characteristic (infectious and immune) 
times. The formulation is made in terms of delay equation which are both numerical 
integrated and linearized. Results from simulation are included showing where they 
support the linear analysis and explaining why not where they do not. Considerations 
and comparison with real diseases are presented along. 

PACS numbers: 87.19.X-, 87.23. Cc 
1. Introduction 

Many diseases that have affected and still affect humans come and go with time in a well 
established way. Examples are plenty and fill the bulletins of world and national health 
organizations. Measles, typhus and cholera epidemic waves, just to cite a few, are even 
part of mathematical biology books [U [2] . The common denominator of such diseases is 
the cyclic natural history of them, in which a susceptible subject can go to infected, then 
to removed, and finally back to the susceptible state. However, the mere cyclic nature of 
the disease does not grant an oscillatory behavior of its epidemic, as can be exemplified 
by gonorrhea 0,11]. In 2009 a new variant of influenza A H1N1, dubbed swine flu, 
appeared in the scene taking the media to discuss on the wave behavior of influenza. 
In which way do these oscillations arise in a population, apparently synchronizing the 
infective state of many individuals? Are they related to external driving causes, such as 
the seasons? Or do they arise dynamically from the very natural history of the disease? 
It is known that several causes can produce oscillations in model epidemic systems: 
seasonal driving [5], stochastic dynamics [HIE], a complex network of contacts [BJ, etc. 

In every infectious disease several characteristic times clearly appear in the 
dynamics. In a SIRS type disease (see figure [T]) one has an infectious time (during 
which the agent stays in category /, infected and infectious), and an immune time 
(during which the agent stays in category R, recovered from infection, and immune to 
re-infection until it returns to the susceptible class S) . A standard SIRS model uses the 
inverse of these characteristic times as rates in mass-action equations, showing damped 
oscillations toward the endemic state in most typical situations. However, the standard 
SIRS model does not exhibit sustained oscillations for any value of the parameters. 
Would it be possible for a deterministic SIRS model to sustain oscillations, and if so 
how do the characteristic times relate to the period of the epidemic? These are some of 
the questions we address in this contribution. 
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As observed by Anderson and May [T] the mathematically convenient treatment of 
the duration of the infection as a constant rate is rarely realistic. It is more common 
that recovery from infection takes place after some rather well defined time. It would 
seem a valid simplification to assume that recovery happens exactly at a given time, 
instead of continuously at some rate. Both extremes are called Type B and Type A 
recovery respectively by Anderson and May. Most infections belong to an intermediate 
type between these two extremes (but closer to Type B). These intermediate types have 
an infectious time distributed with a shape between that of an exponential (Type A) 
and that of a delta distribution (Type B). An accurate model should implement the 
distributions based on empirical data. And as suggested by Hoppenstaedt and others 
[9], that could be done replacing the simple constant rate term by an integral, leading to 
integro-differential equations. Similar considerations can be made about the transition 
from recovered to susceptible term. It is almost thirty years ago that Hethcote et al [TU] 
showed precisely that the SIRS model with fixed delay in the recovered to susceptible 
transition presents stable periodic solution for certain parameter values. On the other 
side they demonstrated that fixed delay in the recovery from infection term does not 
have the same effect when the other one is considered in the usual way. This implies 
that the immunity time plays the most important role in the emerging oscillations as we 
will show in detail. Recently, continuing with the work of Hethcote et al , Taylor and 
Carr [TT] studied in detail the dynamics of the SIRS model with temporary immunity, 
but considering that a fraction of the population acquire permanent immunity. In other 
words, that is equivalent to a distribution of immunity times made of two delta peaks 
(at a finite time and at infinite). That makes the analysis of the oscillations much 
more involved because of the extra parameter: the fraction of individuals who became 
permanently removed. 

Our contribution can be regarded as an extension of the work of Hethcote et al 
with time delays, in which we use arbitrary distributions of both infectious and immune 
times. While doing this we want to keep the problem as simple as possible in terms 
of parameters, thus we avoid the use of vital dynamics. We start by considering the 
most extreme case: fixed delay in both terms (delta distributions). This is the simplest 
mathematical case in terms of delays, and it is close related with the first case studied 
by Hethcote et al. P2] • Then, we consider a mono-parametric family of models in which 
the times are described by continuous distributions — being possible to go continuously 
from a Type B to a Type A model for example. In other words we go from the SIRS 
model with deterministic delays to the classical constant rates SIRS, taking all the 
intermediate situations in between. In this way, while we cover previous results, we can 
go further considering the most general situation for a SIRS model. 

For delayed models, we analyze the onset of sustained oscillations and characterize 
them with the parameters of the system. Linear analysis, together with numerical 
solutions of the nonlinear model, provide a clear characterization of the phenomenon. 
Stochastic numerical simulations provide further support to our analysis. Moreover, 
we show the effect of the shape of the distributed delays on the stabilization of the 
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oscillations. Besides, the oscillations period, which may play an important role in the 
design of intervention policies, are shown to satisfy general rules in terms of the SIRS 
parameters. 



2. SIRS model with arbitrary recovery and loss of immunity dynamics 

As mentioned above, the usual formulation of a SIRS model implies that recovery from 
infection, / — > R, proceeds at a rate which is independent of the moment of infection. 
Also, the loss of immunity R — > S is also just proportional to the current sub-population 
and proceeds at its own rate. The mathematical formulation of such a model is usually 
presented in terms of differential equations as follows: 

M5) = _^ Wl(i) + ^), (ls) 

at r r 
di M- Mi)i(t )-*), (U) 



dt T, 

dr(t) i{t) r{t) 



:ic) 



dt Ti T r 

where s(t), i(t) and r(t) stand for the corresponding fractions of susceptible, infectious 
and recovered individuals in the population (s(t) + i{t) + r(t) = 1). The parameters 
of the model are /3, the contagion rate per individual, and Tj and r r , the characteristic 
infectious and immunity periods respectivehjj] 

The analysis of more general systems — in which the recovery from infection and 
loss of immunity processes obey more general and more realistic dynamics — is more 
involved, leading to non-local integro-differential equations. Before proceeding to the 
most general situation, we analyze the simplest case of fixed times. 



2.1. SIRS with fixed infectious and immunity times 

Let us assume that the disease is characterized by an infectious time Ti as well as 
an immunity time r r . That is, an individual that becomes infectious at time t will 
deterministically recover at time t + Tj, becoming immune, and will subsequently loose 
its immunity at time t + Tj + r r = t + r , becoming susceptible again. The process is 
schematically depicted in figure [TJ This system can be represented by the following set 

SIRS 
1 1 1 * 

/« Xj wn X R M time 

contagion 

Figure 1. Timeline of an individual, showing the course of the disease after contagion. 

| Note on nomenclature: the infectious time is frequently called as recovery time as well, because it 
marks the passage from I to R which is the recovery from infection. While we prefer the first name, 
in order to avoid ambiguities with labels we associate them to the classes, i.e.: is the time in the 
infective class, r r is the time in the recovered class. 
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of equations for the fraction of susceptible and infectious sub-populations (bear in mind 
that r(t) = 1 — s(t) — i(t), so that just two equations describe the dynamics): 
ds(t) 



dt 
di(t) 

~dT 



■ -ps(t)i(t) + ps(t-r )i(t-ro), 
Ps{t)i{t)-Ps{t-Ti)i{t-T^. 



(2a) 
(2b) 



Before proceeding with the detailed analysis of the above formulation of the SIRS 
model we show in figure [2] a result in advance, comparing the time evolution of the 
fraction of infectives in the two scenarios: the SIRS with two fixed time delays and the 
standard SIRS. The last one produce the well know behavior of damped oscillations 
toward the endemic state. Using the same parameters (/3 = 0.4, Tj = 5, r r = 50) the 
SIRS with two delays shows clearly the sustained peaked oscillation in the infective 
fraction. The time delay in the removed to susceptible transition instead of the usual 
continuous rate transition is the responsible for the oscillation as we will see in the 
next sections. In general, delay equations applied to SIR like models have noticeable 
effects on the dynamics, that is different forms of the infection time distribution yield 
different dynamics for the infectives and susceptibles [T2]. The example of figure [2] is a 
remarkable case. 
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Figure 2. Time evolution of the fraction of infected individuals for the SIRS model. 
Comparison between the numerical solutions for the standard and the two fixed time 
delays formulation. Parameters are j3 = 0.4, Tj = 5,r r = 50. 



In equations (I2all2 b\\ . the first terms represents the contagion of susceptible by 
infectious ones, which occurs locally in time at a rate /3. The second terms account 
for loss of infectivity ( |2fej) and loss of immunity ( 1211 . Both terms correspond to the 
individuals infected at some earlier time: t — and t — tq respectively, and who 
have proceeded through the corresponding stage of the disease. These equations must 
be supplemented with initial conditions, appropriate for interesting epidemiological 
situations. A reasonable choice, which we use in the remaining of the paper, is an 
introduction of infectious subjects into a completely susceptible population: 

s(0) = 1 - zo, z(0) = io, r(0) = 0. (3) 
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The system ( 12 atl2B has the drawback that any pair of constants (so, ^o) satisfies 
them, apparently indicating that any pair of values are equilibria. The origin of this 
problem lies in the fact that ( l2aH26j) together with (jHJ) do not constitute a well-posed 
differential problem. Due to the non-locality in time, extended initial conditions must 
be provided. Mathematically, it is usual to provide arbitrary functions s(t) and i{t) 
in the interval [— Tq,0). From an epidemiological point of view, however, it is more 
reasonable to provide just the initial conditions at t — 0, and complementary dynamics 
in the intervals [0, 7*): no loss of infectivity or immunity, just local contagion; and [rj, To): 
transitions from I to R (the second term of (12 ap being absent), and the functions s(t),i(t) 
already obtained by the initial dynamics. 

Indeed, this is the most reasonable choice for the numerical solution of the system, 
and it is the one we have followed in the numerical results shown below. For the analysis 
of equilibria, however, an integral representation of the system results into a better-posed 
problem and the difficulties for the calculation of the equilibria disappear. 

An integral equation equivalent to ( j2jjj) is: 



the interpretation of which is immediate: the integral sums over all the individuals that 
got infected since time t — Ti up to time t. These are all the infectious at time t, since 
those infected before have already recovered. The integration constant c\ is, in principle, 
arbitrary, but it is easy to see that it must be zero since no other sources of infectious 
exist beyond those taken into account by the integral term. 

Complementing (jlj) it is convenient to write the equation for 1 — r = s + i, which 
cancels out the first terms in (12 aH2 b\i : 



where, again, C2 is an integration constant. In this case we have C2 = 1 since no other 
sources of R exist. 

The system ( )4ll5i) can be solved for the equilibria of the dynamics, s* and i*. One 
obtains: 



which coincides with the equilibria found numerically by integrating ( l2~aH2~B . and also 
corresponds to the same equilibria that can be found in a constant-rate SIRS model 
([Taj) , where the rates of recovery and loss of immunity are 1/n and l/r r respectively. 

2.2. SIRS with general distribution of infectious and immunity times 

The idea behind the fixed-time delays can be generalized to describe more complex 
dynamics. Let us start with the infected individuals, which is simpler. Consider a 
probability distribution function G(t), representing the probability (per unit time) of 
loosing infectivity at time t after having become infected at time 0. Observe that the 




(4) 




(5) 
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fixed-time dynamics is included in this description, when G(t) = 5(t — Tj). G(t) can 
be used as an integration kernel in a delayed equation for the infectious. Indeed, the 
individuals that got infected at any time u < t and cease to be infectious at time t are: 



(3 / s(u)i(u) G(t — u)du, 
Jo 

so that the differential- delayed equation for i(t) is: 

^ = ps(t)i(t) -/3 J* s{u)i{u) G{t - u)du, (7) 
analogous to fl2fe|) . 

A second kernel H(t) must be considered for the loss of immunity process. The 
differential equation for susceptible can then be written as: 
ds(t) i-t r rv 



dt 



-ps(t)i(t) + /3 



s(u)i(u) G{v — u)du 



H{t - v)dv, 



where the second term corresponds to individuals that get infected at earlier times, then 
loose their infectivity at intermediate times with probability G, and finally recover at 
time t with probability H. 

The difficulty with initial conditions that we faced in the fixed-times system is also 
found here, and can be solved in the same way. Integral equations for i(t) and s(t) +i(t) 
result: 



i(t) — ci + (3 / s(u)i(u)du — (3 



s(u)i(u) G(v — u)du 



dv, 



s(t)+i(t) = c 2 + (3 



JO 



s(u)i(u)G(v — u)duH(x — v) 



—s(v)i(v)G(x — v) 



dv dx. 



(9) 



(10) 



These equations can be used to find closed expressions for the equilibria which, 
depending on the functional forms of G and H, can be solved analytically. In general 
one finds two sets of solutions, the disease free one: s* — 1, i* — 0, and the endemic one 
bifurcating from it: 

1 .* 



/3(S!-E 2 ) : 



with: 



Ei = / 1 - / G(v - u)du 
Jo L Jo 



dv, 



o Jo 



H(x — v) j G(v — u)du — G(x — v) 



dv dx. 



(11) 

(12) 
(13) 



We observe that using either Dirac deltas or exponential functions for the kernels 
(corresponding to fixed-times and constant rates, respectively), these integrals can be 
solved analytically to find the equilibria. 

Between the two extremes of constant rates and fixed delays, as mentioned, realistic 
systems are expected to display some extended probabilities distributions for infectious 
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and immunity times. A convenient interpolation between the exponential and delta 
distributions that characterize those regimes can be achieved by gamma distributions: 

Pi fPi-lp-pit/n 

G » (t) = V'te-D! ■ (14a) 
H ^ = % P ,.-iy. ■ < 146) 

These distributions have mean Tj and r r respectively, for any value of the parameter 
Pi ir . Besides, they interpolate between exponential (when p^ r = 1) and Dirac delta 
distributions (when p i>r — > oo), with smooth bell-shaped functions for intermediate 
values of p i>r . It can be shown that the equilibria (TTTT) are identical to the classical ones 
(EJ) for any p i>r . 

3. Sustained oscillations in SIRS with delays 

Standard SIRS systems, without delays (or equivalently, with Gi(t) and Hi(t) as delay 
kernels) have either nodes or stable spirals as equilibria. That is, oscillations appear in 
them as transient regimes damped towards the fixed points. SIRS systems with delays, 
on the other hand, can exhibit sustained oscillations. These appear as a Hopf bifurcation 
of the spiral points, controlled by the parameters Tj, r r and (3. A linear stability analysis 
of the fixed-times case can exemplify how this happens. 

Assuming that the system (l2aH26l) is close to equilibrium, one sets s(t) = s* +x(t), 
i(t) = i* + y{t), obtaining in linear approximation a linear delay-differential system for 
the departures from equilibrium: 

x(t)/(3 = -i*x(t) - s*y{t) + i*x(t - r ) + s*y{t - r ), (15a) 
y(t)/P = i*x(t) + s*y(t) - i*x(t - n) - s*y(t - n). (155) 

From this system, proposing exponential solutions x(t) = c\e xt and y(t) = C2e xt , a 
transcendental characteristic equation is obtained: 

A 2 + A/3 [s*(e~ Xn - 1) - i*(e' Xro - 1)] = 0. (16) 

Equation ( Tl6|) can be solved numerically for complex A, obtaining from its real part 
the bifurcation line from the stable spirals. This line is shown in figure [3] along with 
the amplitude of oscillations represented by a colour (gray) map. The amplitude is 
the result of the numerical integration of the full nonlinear system (I2all26|) . Figure [3] 
condenses the bifurcation phenomenon as a function of the key parameters r r /ri and 
Rq. The black region represents the non oscillating endemic solution. It can be seen 
that the linear analysis, represented by the line, defines almost exactly the transition. 

3.1. Distributed delays: the general case 

The fact that we have qualitatively different results regarding the nature of the endemic 
state for a delta or an exponential distribution gives rise to a fundamental question. Is 
the existence of an oscillatory endemic state particular to the delta distribution? Or is 
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0.5 1 1.5 2 2.5 3 3.5 
T r /Ti 



Figure 3. Bifurcation diagram of the SIRS model with fixed times in the space 
defined by T r /ri and R . The ycllow(white) line shows the linear result. The squared 
root amplitude of the infectives oscillations is shown as colour (gray) coded shades 
above the transition line. The black region means zero amplitude, representing non- 
oscillatory endemic states. 



there a critical shape of the distributions G and H necessary for the emergence of the 
oscillations? Using the Gamma functions defined in ( \14ah we can check the existence 
of such solutions for different shapes by controlling the parameter p i>r . In this way, we 
can verify if there is a critical shape pi^ r = p c ((3,Ti,r r ) beyond which the system has 
sustained oscillations at the endemic state. Linearizing the general system (171151) in the 
same way presented in the previous section, one has the following integral characteristic 
equation: 



A 2 + \pe 



o Jo 



t-v 



H(v)G(u)e- x{u+v) dudv 



A/3s* 



G(u)e~ du 



0. 



;i7) 



Using the distributions fll4ap in the equation above and taking the limit t — > oo we get: 



A 2 + A/3f 





-Pr 



- A/3s* 




0. 



As expected, for p itr = 1 and pi^ r — > oo we recover the characteristic equations of 
the constant rates and fixed delays models, respectively. Solving (1181) numerically for 
complex A, one finds that for every set of parameters (provided that Ro = f$Ti > 1) 
where p r > 1 there is always a critical shape pi = p c above which the endemic state 
consists of sustained oscillations. Conversely, any value of p r > 1 can present sustained 
oscillations in some region of parameter space. Instead, the constant rate model p r = 1 
is a particular case where there is no such solution for any set of parameters, as was 
demonstrated by Hethcote et al ([TU]). 

Figure H] shows the real part of the solution of (fT8j) as a function of the shape 
(Pi,r — p) of the two time distributions, and for fixed values of the parameters f3, and 
T r . In other words we can appreciate (for those j3, r i: r r ) the critical p value (p c ) to have 
oscillations in the SIRS dynamics, i.e. the value at which Re[\] = 0. 
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Figure 4. Real part of the eigenvalue A as a function of the shape parameter 
Pr — Pi — P- Inset: examples of the distribution G p with tj = 1. 



The critical shape p c obtained by numerical calculation of ( Tl8|) gives an accurate 
prediction of the emergence of oscillations in the full nonlinear system (ITfSj) . as can be 
seen in figure EJ The meaning of the critical shape is simple: the time distributions have 
to be narrower than the one represented by p c in order to have sustained oscillations. 




2000 2020 2040 2060 2080 2100 

time (arbitrary units) 



Figure 5. Numerical integration of the SIRS general model with j3 = 1, t; = 5, 7> = 10 
near the critical shape pi_ r = p c predicted by the linear analysis of the system (inset, 
12 <p c < 13). 



The distributions of infectious and immunity times used so far share a common 
shape given by the value of p — that made the analysis, restricted to only one shape 
parameter, simpler. In real situations, though, it is reasonable to expect that these 
uncorrelated kernels have different shapes, not necessarily of the same relative width. 
We explore this more general scenario, presenting an oscillation diagram in terms of pi 
and p r in figure El for two sets of the SIRS parameters. Some interesting conclusion can 
be extracted from such diagram: 

• If the immunity time distribution is not sharp (p r < 5 for the parameters used in 
figure E]) there is no oscillation independent of the type of infection time distribution. 
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• For a relatively narrow distribution of immunity times there is a critical infective 
time width above which (i.e. a critical Pi below which) no oscillations are obtained. 

• The broader one of the time distribution is, the narrower the other one has to be 
in order to have oscillations. 

• Longer immunity times (in units of infectious time) makes the oscillatory region 
wider in the Pi,p r space. 
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Figure 6. Oscillation diagram in the pi, p r plane for two sets of parameters: /3 = 1, 
Tj = 2, r r = 10 (circles) and /3 = 1, Tj = 2, r r = 20 (triangles). 



In the case Pi yT — > oo there is a critical value of Ro(T r /Ti) above which the endemic 
solution is always a stable cycle. The linear analysis also shows that for finite pi >r the 
bifurcation is more involved. Figure [7J shows that the bifurcation line encloses a region 
of oscillating solutions. Then, for a given value of r r /rj, there is a second critical Ro, 
larger then the previous one, where the endemic solution ceases to be cyclic. This 
phenomenon is verified in the numerical solution of the nonlinear system. For any given 
r r /Ti, if one increases p the region of oscillation grows (as shown in figure Ej), so that 
in the limit p i r — y oo the upper critical Ro disappears. In such a situation, the only 
way to break the oscillations is by decreasing Rq. On the other hand, by decreasing 
Pi, r the oscillatory region shrinks, disappearing completely when p iyT = 1 (exponential 
distributions, constant rates). 

Yet, the interplay between the uncorrelated shapes pi and p r and the SIRS 
parameters is very interesting. Again, in figure [7J we see for example that when 
Pi = Pr = 5 (equal very broad shapes) the oscillation region is very small. But increasing 
Pi to 20 (narrow infective times distribution), keeping the other fixed, expands the region 
considerably. Similarly, but in other part of the same diagram, when pi = p r = 20 
(equal very narrow shapes), the region of sustained oscillations is relatively large. Now, 
by letting the distribution of infective times to be broader (p^ = 5) again with the 
other distribution fixed, we see the oscillatory region to shrink significatively. So, while 
it is true that the immunity times distribution must be necessarily different from an 
exponential in order to have oscillations, in general both distributions are relevant to 
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define the actual parameter region of oscillations. The narrower one or the two are, the 
wider the oscillatory region is. 

10 



8 
6 
4 



2 



"0 5 10 15 20 

X I X. 

r i 

Figure 7. Oscillation diagram in the Tv/ri, Rq plane for the SIRS model with different 
distribution for infectious and immunity times controlled by the p;. r -shapc factor of the 
Gamma distribution function. For each pair, pi, p r the oscillatory region is enclosed 
by the corresponding line. The curves with only one p are for Pi = p r = P- 

Three examples of dynamics of the SIRS with different shapes sets are displayed 
in figure [SJ The curves correspond to the numerical solution of systems with the same 
epidemic parameters but three different shape sets: (pi = p r = oo), (pi = l,p r = oo), 
and (pi = oo, p r = 1). The first two show the sustained oscillations because the SIRS 
parameters are inside the oscillatory region, while in the last case we know that there is 
no possible region of oscillation so we have the damped behavior. Besides, we see that 
enlarging the Tj distribution alone, while it shortens the period of the oscillation, it does 
not destabilizes the oscillation. Remarkably it shows an enhancement of the number of 
infected during the low part of the cycles. This behavior diminishes the probability of 
extinction in a discrete population realization of the system, a problem that is common 
in simulations as we will see in the next section. 

From the imaginary part of (jl8p it is possible to obtain the period of oscillations in 
the linear approximation. In figure [9] we show the result of this calculation for different 
shapes of the distributions of infectious and immunity times. For finite p^ r , the period 
of oscillation has a dependency on the parameters that approaches the one found for 
Pi X — » oo (delta distributions, fixed times) in the lower critical value of Rq, given by 

T = 3/4r; + 2r r . (19) 

For a fixed value of r r /rj, this period decreases as Rq increases up to the upper limit 
of oscillation, being bounded from below by 3/4rj + r r , as shown in figure There 
are no oscillations above and below these two straight lines, for any value of pi sr . This 
remarkable lock of the period of disease oscillation is related to the one obtained by 
Taylor and Carr ([TT]) in the case of exponential distributed studied by them. 
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Figure 8. Fraction of infectious as a function of time for three distributions of r r and 
Tj with different shapes but same parameters (3 = 0.4, Tj = 5, r r = 50. 
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Figure 9. Period of oscillations as a function of r r /r,;. Each curve corresponds to a 
different shape of the time distributions, as shown in the legend where Pi = p r = P- 
There are no oscillations above and below the dashed lines. 



4. Simulation 



A complete picture of the general SIRS dynamics needs to contemplate straightforward 
numerical simulations. We believe that simulations represent the most accurate 
implementation of the real system, which is discrete and stochastic. Therefore, as a test 
of our analytical and numerical results with the generalized SIRS model, we present 
here results from a probabilistic discrete model. In this model a finite number of agents 
meet at random and contagion proceeds in a probabilistic way. The phase diagram for 
such a system, with delta distributions for r» and r r , is shown in figure (TDJ The temporal 
evolution of the system is shown in figure [TT] along with the numerical solution of the 
deterministic model. 

The agreement between the two implementations of the same model is very 
satisfactory, as can be seen by both the bifurcation diagram and the temporal solution. 
Regarding the diagram the agreement is not as perfect as with the numerical solution 
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Figure 10. Diagram of the oscillation amplitude, in the T r /Tj, Rq plane, obtained by 
simulation with N = 1000 agents in the SIRS model with fixed times. Vertical axis is 
Rq = (3ri, horizontal axis is the ratio r r /Tj, and the colour (gray) map represent the 
squared root amplitude of the i(t) at the steady state. The black region means zero 
amplitude, representing non-oscillatory endemic or non-epidemic states. Superimposed 
is the critical line obtained by the linear analysis as in figure 
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Figure 11. Evolution of an epidemic with stochastic and deterministic models. Time 
distributions are deltas at Tj = lOdays and r r = I80days, and ft = .ISperday. With 
the time scale in years, the annual oscillation is clearly observed (T = 2r r + 3/4ri = 
367days). 



of the deterministic equations (figure [3]) because of fluctuations that arise in simulation 
due to finite size. These cause large amplitude oscillations which become extinct but 
that contribute with a non-zero value for the amplitude. As for the temporal evolution, 
eventually, as time goes by, a small drift develops in the simulation due to the stochastic 
nature of its dynamics. Far from being disappointing this is a desirable feature, since in 
real systems epidemic oscillations are not exactly periodic. Yet, the extinctions, that are 
normally observed in the simulations, are related with the discreteness of them. Long 
time ago discussed by Bartlett [13] and others it manifests here as being very sensitive 
to the distance to the lower threshold in the bifurcation diagram, which relates to the 
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amplitude of the oscillation. Greater amplitudes drive the system closer to the absorbing 
state at i — 0. To illustrate this effect we show in figure [15] three dynamics obtained by 
simulation with N = 10 5 for three different values of (3, close to the onset of epidemics 
and oscillations. However we saw in the previous section (figure [H]) that a distribution 
of infectious times with finite width can increase the lower values of the infectives cycle, 
which eventually may prevent the extinction in the simulation counterpart. 
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Figure 12. Evolution of the epidemic in the stochastic model, showing extinctions of 
the infected population. Time distributions arc deltas at n — 10 and r r = 180 days, 
N = 10 5 . 




5. Conclusions 

We have analyzed a general SIRS epidemic model, in which the infective state as well 
as the immune state last prescribed times drawn from distributions. These distributed- 
time transitions constitute a generalization of the standard SIRS model, in which the 
transitions out of the infective and the immune states happen at a constant rate. 
The generalization allows for situations in which these states last for certain fixed 
times — which is more similar to many real diseases than the constant rate assumption. 
Between the two extremes of constant rates and fixed times, we have also analyzed the 
intermediate situations of broader or narrower distributions of the transition times. 

The generalization runs along the proposals made by previous authors, by 
implementing the differential equations of the model as non-local in time. For example 
Hethcote etal. have shown that cyclic models have a transition to oscillatory behavior 
when the immunity time is distributed [TO] . 

Our contribution shows how the oscillating state arises as a function of all the 
parameters of the model, completing a phase diagram that provides a thorough and 
general view of the possible behaviors. We show, moreover, that the linear analysis, the 
numerical integration of the model, and its stochastic implementation, all converge to 
the same general picture, within the inherent limitations of each one. 
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For delta distributed delays (i.e., for fixed times spent in the infective and recovered 
classes), the phase diagram of figure [3] shows that the region of oscillation is bounded 
from below by the basic reproductive number Rq as a function of the scaled immunity 
time T r /rj. We see that, as r r decreases toward zero, the minimum value of Rq diverges. 
This result is in agreement with the fact that in the SIS model there are no sustained 
oscillations. On the other hand, the amplitude of the oscillations grows both with Rq 
and T r /Tj. This kind of behavior may be relevant in the analysis of real epidemics where 
changes in the parameters are occurring due to interventions, advances in treatment or 
natural causes. For example, let us imagine an epidemic in the endemic non-oscillating 
region (shaded black in the diagram), for which, as a desirable consequence of the 
treatment of the disease, there is an increase in the duration of the immune state, 
at constant Rq. As a result, the epidemic may start to oscillate. Such a transition 
may manifest itself as an initial increase of the infectious fraction, due to the onset of 
oscillation, which should be properly understood in the right context. 

We have also shown that the transition to oscillations behaves as a critical 
phenomenon depending on the width of the time distributions. Figure H] shows 
this behavior. The remarkable stabilization of the oscillations of pertussis after the 
introduction of massive vaccination [TJ section 6.4.2] may be related to a narrowing of 
the immunity time distributions. 

For distributed delays, corresponding in our model to values of the parameter 
1 < Pi,r < oo, we have found a remarkable reentrance phenomenon in the phase diagram. 
As exemplified by figure [7J this reentrance means that the region of oscillations is also 
bounded from above by a curve of Rq vs T r/ T i- One sees, also, that the region of 
oscillation shrinks with lowering p^ r , until its disappearance at the critical value p c . 
Only the case p i>r — > oo, corresponding to delta distributed times is unbounded. This 
is at variance with the SIR case analyzed by [14] who claims that there is no significant 
change between large values of p^ r = p (p = 10-20) and p — > oo. On the other hand our 
model supports their result in that the region of oscillations is already large for p in this 
range. An important general result that can be derived from this diagram is the fact 
that, whenever r r ^> Ti, the minimum value of Rq is very close to 1, implying that such 
systems — provided that pi )T is sufficiently large — are very prone to sustained oscillations. 
This implies that oscillations due to the implementation of real time distribution into 
the SIRS formalism are not a marginal or delicate effect as someone could think before. 

Our analysis also shows the behavior of the period of oscillation and its dependence 
on system parameters. Is interesting to discuss this aspect of our model in connection 
with real infectious diseases. For example in the case of influenza, where the period of 
oscillation is more or less a year and the infectious time is of the order of 10 days, our 
model predicts oscillations for any Rq > -Romm = 1-09. This minimum value is lower 
than the usual estimate for this disease, which is around 1.4 [15]. Indeed, this value is 
so close to Rq — 1 that one can conclude that there will always be oscillations for this 
disease as long as it remains endemic. Besides we see that using a value of r r = 180 
days for the loss of immunity, the region of oscillation lays between 1.22 < R < 2.4 
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for a rather wide distribution with pi )T — p = 10, and 1.16 < Rq < 3.4 for a narrower 
distribution with p = 30. In the first case, the corresponding period of oscillation is 
188 < T < 342 (in days), while on the second one it is 190 < T < 360. Therefore, one 
needs a narrow distribution to satisfy the observed period of oscillation for influenza. In 
any case, the predicted region of oscillation includes the observed values of Rq. We want 
to stress that, although an infection by a specific strain of the influenza virus confers 
permanent immunity after the infection period, the mutation rate of the virus can be 
thought of as an effective immunity period and therefore can also be modeled with the 
SIRS dynamics. 

In the case of other oscillating epidemics, such as syphilis and pertussis, our model 
also fits the observed data rather well. For syphilis, in the p^ r — >• oo case, using the 
period of 132 months measured in United States [3] and the usual value of = 6 months, 
we have an Ro m in = 1-17 and r r « 6 years. Greater values of Rq (which are expected 
for syphilis in [3]), with fixed Tj, predict even longer immune times (which is also the 
case of syphilis). For pertussis, using the period of 4 years and Tj = 1 — 2 months, we 
find Ro m i n = 1.06 — 1.15, which is smaller than the estimations obtained for mixing 
models [T5] . 

In the open field of oscillatory epidemics, we believe that there is a multiplicity of 
causes that can concur to produce the observed phenomena. Our present contribution 
does not intend to provide an exclusive and definitive answer to this matter. The 
construction of valid theoretical models for real diseases should incorporate all the 
relevant mechanisms, therefore a thorough theoretical investigation of any concurrent 
possible cause deserves due attention. We hope our work has revealed one good 
candidate. 
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